t=550; % t is the period from which one starts the projection
max_rep=500; % number of times calculations are repeated for a different row of std res eps
muu=load('mu_all3.txt');
eps=load('eps_all3.txt');
hilos=load('hilos.txt');
sper=length(eps(:,1));
eps_s=sort(eps);
giorni1=load('date_7_18_95.txt');
MATLABDate = x2mdate(giorni1);
giorni=datestr(MATLABDate,2);
% load matrices
aplus=load('aplus.txt');   % impact matrix lag 1 overall period dpos
aminus=load('aminus.txt'); % impact matrix lag 1 overall period dneg
acplus=load('acplus.txt');   % impact matrix lag 1 crisis period dpos
acminus=load('acminus.txt'); % impact matrix lag 1 crisis period dneg
ac2=load('a_lag2.txt'); % impact matrix lag 2 
omega1=load('om.txt');
dneg=load('dneg.txt');
beta=load('beta.txt');
ac1=beta+(aplus+aminus)/2;
ad1=(acplus+acminus)/2;

% from here down the same as before

c_e=cov(eps); % covariance matrix of eps
sd=sqrt(diag(c_e));
k1=200; % total length
k=100; % for display purposes
% select a row from muu
if ge(t,473) & le(t,842) % crisis period
   ac=ac1+ad1;
   om=omega1(:,1)+omega1(:,2);
elseif ge(t,843) % post crisis period
   ac=ac1;
   om=omega1(:,1)+omega1(:,3);
else
   ac=ac1;
   om=omega1(:,1);
end
mu=muu(t,:);
mu=mu';
nm=length(mu);
mu=[mu ; hilos(t-1,:)'];
a=[ac ac2; eye(nm) zeros(nm,nm)];
om=[om ; zeros(nm,1)];

% end of initializations

halflife=zeros(nm,nm);
peaks=zeros(nm,nm);
troughs=zeros(nm,nm);
max_2=zeros(nm,nm);
b=zeros(nm*2,k1,nm);
b1=zeros(nm*2,k1,nm); % baseline 
b2=zeros(nm*2,k1,nm); % shocked solution
b3=zeros(nm*2,k1,nm); % ratio
s_mat=zeros(nm,nm);
for mkt=1:nm
    if mkt ==1
        s=[sd(mkt) ; c_e(2:nm,mkt)/sd(mkt)];
    elseif mkt==nm
        s=[c_e(1:nm-1,mkt)/sd(mkt) ; sd(mkt)];
    else
        s=[c_e(1:mkt-1,mkt)/sd(mkt); sd(mkt); c_e(mkt+1:nm,mkt)/sd(mkt)];
    end 
    s_mat(:,mkt)=s;
    s=[s;zeros(nm,1)];
    b(:,1,mkt)=a*(s.*mu); % the vector is now multiplied by mu element by element
    b1(:,1,mkt)=om+a*mu;
    b2(:,1,mkt)=om+a*((1+s).*mu); 
    for i=2:k1
         b(:,i,mkt)=a*b(:,i-1,mkt);
         b1(:,i,mkt)=om+a*b1(:,i-1,mkt);
         b2(:,i,mkt)=om+a*b2(:,i-1,mkt);
    end
end
b=b*log(10)*sqrt(252)*sqrt(pi)/sqrt(8);
b1=b1*log(10)*sqrt(252)*sqrt(pi)/sqrt(8);
b2=b2*log(10)*sqrt(252)*sqrt(pi)/sqrt(8);
b3=b2./b1;


%%START GRAPHS - input b matrix which is nm,k1,nm first dimension contains from by
%%row last dimension contains "to" by row
tt=(1:k);
ctry=['HK';'IN';'KO';'MA';'PH';'SI';'TA';'TH'];
mkt_mrr=['h';'o';'*';'^';'+';'p';'s';'p'];
%mkt_mrr=['none';'none';'none';'none';'none';'none';'none';'none'];
mkt_col=['k';'r';'g';'b';'k';'m';'y';'c'];
mkt_col=['k';'r';'g';'b';'k';'m';'b';'r'];
mkt_sty=['-   ';'-   ';'--  ';':   ';'-.  ';'-   ';'-   ';'-   '];
mf = figure(1);
set(mf,'Color','w')
%annotation('textbox',[0.0 0.85 0.7 0.1],'String',['Conditional Response to a 1sd shock in a market. Starting Date: ' giorni(t,:)],'FontSize',16) 
%title(['Conditional Response to a 1sd shock in a market. Starting Date: ' giorni(t,:)],'FontSize',16 )
%coo=[0.05 0.5 0.4 0.35 ; 0.55 0.5 0.4 0.35 ; 0.05 0.05 0.4 0.35 ; 0.55 0.05 0.4 0.35]; 
plot(tt,b3(1,1:k,1)-1,'LineWidth',2,'MarkerSize',6,'Marker',mkt_mrr(1,:),'Color',mkt_col(1,:));
title(['Originating Market: ' ctry(1,:) ' - Starting Date: ' giorni(t,:)] ,'FontSize',24)
hold on
for ii=2:8
   plot(tt,b3(ii,1:k,1)-1,'LineWidth',2,'MarkerSize',6,'Marker',mkt_mrr(ii,:),'Color',mkt_col(ii,:));
end
hold off
legend(ctry,'Location','NorthEast')
xlabel('Horizon','FontSize',24)
ylabel('Response','FontSize',24)
ylim([0 .45])
xlim([1 100])
set(gca,'FontSize',24)
X=getframe(gcf);
imwrite(X.cdata,[ctry(1,:) '_IRFOct22_DD.png'],'png','BitDepth',16);
